//cubic spline algorithm is
//based on 131page of AOKI's thesis 

#include "gnuplot_i.hpp"
#include <iostream>
#include <fstream>
#include <stdlib.h>

double newton(double x[], double y[], int n ,double t){
  static int flag = 1;
  static double a[100];
  double w[100], s;

  int i,j;

  if (flag == 1){
    for (i=0; i<n; i++){
      w[i] = y[i];
      for(j=i-1; j>=0; j--)
	w[j] = (w[j+1]-w[j])/(x[i]-x[j]);
      a[i] = w[0];
    }
    flag=-1;
  }

  s=a[n-1];
  for(i=n-2;i>=0;i--)
    s=s*(t-x[i])+a[i];
  return s;
}


double lagrange(double x[], double y[], int n, double t){
  int i,j;
  double s,p;

  s = 0.0;
  for (i=0 ; i<n; i++){
    p=y[i];
    for(j=0; j<n; j++){
      if (i != j)
	p = p * (t-x[j]) / (x[i]-x[j]); 
    }
    s = s + p;
  }
  return s;
}

//lagrange(x,y,3,t));

double linear_spline(double x[], double y[], int n, double t){
  int i;
  double yt = 0.0;

  for(i=0; i<n; i++)
    if(t < x[0] )
      yt = y[i+1] + ((y[i+1]-y[i]) / (x[i+1]-x[i])) * t;
    else if(t >= x[i] && t < x[i+1])
      yt = y[i+1] + ((y[i+1]-y[i]) / (x[i+1]-x[i])) * t;
    else
      yt = y[i] + ((y[i]-y[i-1]) / (x[i]-x[i-1])) * t;
  return yt;
}

/*
double LinearInterpolate(
   double y1,double y2,
   double mu)
{
   return(y1*(1-mu)+y2*mu);
}
*/

//based "Numerical recipes in C"
void spline(double x[],double y[], int n, double yp1, double ypn, double y2[]){
  int i,k;
  double p, qn, sig, un, *u;

 u=(double *)calloc(sizeof(double), n);


 if(yp1 > 0.99e30)
   y2[0] = u[0] = 0.0;
 else {
   y2[0] = -0.5;
   u[0] = (3.0/(x[1]-x[0]))*((y[1]-y[0])/(x[1]-x[0])-yp1);
 }

 for(i=1;i<n-1;i++){
   sig = (x[i]-x[i-1])/(x[i+1]-x[i-1]);
   p = sig*y2[i-1]+2.0;
   y2[i] = (sig-1.0)/p;
   u[i] = (y[i+1]-y[i])/(x[i+1]-x[i])-(y[i]-y[i-1])/(x[i]-x[i-1]);
   u[i] = (6.0*u[i]/(x[i+1]-x[i-1])-sig*u[i-1])/p;
 }

 if(ypn>0.99e30)
   qn = un = 0.0;
 else{
   qn = 0.5;
   un = (3.0/(x[n-1]-x[n-2]))*(ypn-(y[n-1]-y[n-2])/(x[n-1]-x[n-2]));
 }

 y2[n-1] = (un-qn*u[n-2])/(qn*y2[n-2]+1.0);
 for(k=n-2;k>1;k--)
   y2[k] = y2[k] * y2[k+1] + u[k];

 free(u);
}


int main(){

  Gnuplot gp("points");

  double fa[4];
  double a, b, c, d;
  char equation[100];

  double *x, *y, *y2;
  int n, i;
  double t;

  n = 3;

  x = new double[n];
  y = new double[n];
  y2 = new double[n];

  x[0] = 0;   x[1] = 1;   x[2] = 2;
  y[0] = 0;   y[1] = 2;   y[2] = 1.5;
  
  /*
  for(t=0.0;t<=7.0; t=t+.5)
    printf("%7.2f %7.2f\n", t, lagrange(x,y,3,t));

  for(t=0.0;t<=7.0; t=t+.5)
    printf("%7.2f %7.2f\n", t, newton(x,y,3,t));
  */
  for(t=0.0;t<=7.0; t=t+.5)
    printf("%7.2f %7.2f\n", t, linear_spline(x, y, 3, t));

      //  spline(x, y, 3, 1, 1, y2);

  //  gp.reset_all();

  //  gp.plot_xy(x, y, "points");

  //  gp.set_style("lines");
  //  gp.set_xrange(x[0],x[1]).set_samples(40);
  //  gp.plot_equation(equation, "cubic");

  getchar();

  delete [] x;
  delete [] y;
  delete [] y2;

  return 0;
}


  
  
//  sprintf(equation,"%f*x*x*x + %f*x*x + %f*x +%f", a, b, c, d);
//  printf("%s", equation);


  
